Setup

This script uses the flowmix package from github repository commit number 251b669443e77c30a4519ce25add277fb46718d4. (https://github.com/sangwon-hyun/flowmix/commit/251b669443e77c30a4519ce25add277fb46718d4)

1d 5-cluster model results

This section contains code to produce the plots for the 1d data 5-cluster analysis in the paper (Figure 7, 12, 13, and Table 2).

First, run the simulations from gradients2-run.R, into a destination directory destin. Then, use this code to summarize the results.

This will save the summary of the results into file.path(destin, "summary.RDS").

In this script, I’m loading this file directly, from ./paper-data/1d-cvres.rds

## Load the summary
## cvres = readRDS(file = file.path(destin, "summary.RDS"))

## Instead, load this summary
cvres_filename = "1d-cvres.rds"
cvres = readRDS(file = file.path(datadir, cvres_filename))

## For back-compatibility; prob was called pie
cvres$bestres$prob = cvres$bestres$pie 
class(cvres$bestres) = "flowmix"

We’ll also load the 1d data from the directory ../paper-data. (Normally, the data is saved to file.path(destin, "meta.Rdata"))

## load(file.path(cvres$destin, "meta.Rdata")) ## This is one way to load the data
dat = readRDS(file = file.path("../paper-data", "MGL1704-hourly-paper-1d-diam.RDS"))
countslist = dat$countslist
ylist = dat$ylist

Figure 12

Figure 12 is produced using this code:

## Plot CV scores
## pdf(file = file.path(figdir, "1d-cvscores.pdf"), width=12, height=5)
plot_cvscore(cvres$cvscore.mat)

## graphics.off()

Figure 7

Figure 7 is produced using this code:

## Plot fitted mean.
## pdf(file = file.path(figdir, "1d-means.pdf"), width=12, height=5)
plot_1d(ylist = ylist, countslist = countslist,
        res = cvres$bestres,
        cex_data = 1) %>%  invisible()
mtext("Diameter", side = 2, line = 2.5, cex = 1.5)

## graphics.off()

Table 2

Table 2 can be produced using this code:

cvres$bestres$alpha %>% t() %>% round(3) %>% as.matrix() %>% knitr::kable() %>% print() ##%>%  xtable::xtable(digits=digits)
clust-1 clust-2 clust-3 clust-4 clust-5
intp -1.345 -1.757 -2.421 -1.854 -2.062
b1 0.000 0.000 0.000 0.000 0.000
b2 0.000 0.000 0.000 0.000 0.000
p1 0.000 0.000 0.000 0.000 0.000
p2 0.000 -0.005 -0.231 0.006 0.000
p3 0.000 0.000 -0.007 0.000 0.000
p4 0.000 0.000 0.000 0.000 0.000
sss 0.000 0.000 0.000 0.000 0.000
sst 0.000 0.000 0.000 0.913 0.000
Fe 0.000 0.000 0.000 0.000 0.000
PP 0.000 0.088 0.000 0.000 0.000
Si 0.000 0.000 0.000 0.000 0.000
NO3 0.000 0.000 0.000 0.000 0.000
CHL 0.000 0.000 0.134 -0.045 0.000
PHYC 0.000 0.000 0.000 0.000 0.000
PO4 0.000 0.000 0.000 0.000 0.000
O2 0.000 0.000 0.000 -0.057 0.000
vgosa 0.000 -0.116 0.000 0.000 0.000
vgos 0.133 0.000 0.000 0.000 0.000
sla 0.000 0.000 0.185 0.000 0.000
ugosa 0.000 0.000 0.000 0.000 0.000
ugos 0.000 0.022 -0.134 -0.086 0.000
wind_stress 0.000 0.000 0.000 0.000 0.011
eastward_wind 0.000 0.000 0.000 0.000 0.000
surface_downward_eastward_stress 0.000 0.000 0.000 0.000 0.000
wind_speed 0.000 0.000 0.000 -0.165 0.000
surface_downward_northward_stress 0.000 0.000 0.000 0.000 0.000
northward_wind 0.138 0.000 0.000 -0.063 0.000
ftle_bw_sla 0.032 -0.094 0.000 0.000 0.000
disp_bw_sla 0.000 0.000 0.000 0.018 0.000
AOU_WOA_clim 0.000 0.000 0.000 0.000 0.000
density_WOA_clim 0.000 0.000 0.000 0.000 0.000
o2sat_WOA_clim 0.000 0.280 0.000 0.000 0.000
oxygen_WOA_clim 0.000 0.000 0.000 0.000 0.000
salinity_WOA_clim 0.000 0.000 0.000 0.000 0.000
conductivity_WOA_clim 0.000 0.000 0.000 0.000 0.000
nitrate_WOA_clim 0.000 -0.636 0.000 0.000 0.000
phosphate_WOA_clim 0.000 0.000 0.158 -0.860 0.000
silicate_WOA_clim -0.015 0.000 0.000 0.000 0.000
par 0.018 0.000 0.000 -0.027 0.000 
cvres$bestres$beta  %>% setNames(paste0("Cluster-", 1:5)) %>% do.call(cbind, .) %>% `colnames<-`(paste0("Cluster-", 1:5)) %>%
  round(3) %>% as.matrix() %>% knitr::kable() %>% print()##%>%  xtable::xtable(digits=digits)
Cluster-1 Cluster-2 Cluster-3 Cluster-4 Cluster-5
intp 0.914 0.082 0.516 -0.326 -0.033
b1 0.014 0.032 0.000 0.046 0.038
b2 0.077 -0.023 0.000 0.000 0.000
p1 0.002 0.000 0.004 0.008 -0.002
p2 0.000 0.006 -0.004 0.013 -0.006
p3 0.007 0.005 0.000 0.010 -0.003
p4 -0.004 0.002 0.010 -0.006 -0.002
sss 0.000 0.000 0.000 -0.067 0.000
sst 0.000 -0.004 -0.028 -0.070 0.000
Fe 0.000 -0.011 0.000 0.001 0.000
PP 0.066 -0.008 0.000 0.000 0.000
Si 0.014 0.000 -0.021 -0.006 0.000
NO3 0.000 0.000 0.000 0.000 0.000
CHL -0.141 0.000 0.000 0.000 0.000
PHYC 0.000 0.000 0.000 0.010 0.000
PO4 0.024 -0.001 0.000 0.000 0.000
O2 0.071 0.000 0.000 0.000 0.000
vgosa 0.000 0.000 0.014 0.000 0.025
vgos 0.012 0.000 0.000 -0.006 0.000
sla 0.004 0.001 -0.012 -0.022 0.002
ugosa 0.000 -0.005 0.000 -0.050 0.007
ugos 0.005 0.000 0.007 0.058 0.000
wind_stress 0.011 -0.029 0.000 0.025 0.000
eastward_wind -0.023 -0.020 0.000 0.021 0.000
surface_downward_eastward_stress 0.000 -0.010 0.000 0.000 0.000
wind_speed 0.000 0.000 0.000 0.000 0.033
surface_downward_northward_stress 0.000 0.000 -0.002 0.011 -0.007
northward_wind 0.000 -0.008 -0.009 0.000 0.000
ftle_bw_sla 0.005 0.024 0.000 -0.011 0.000
disp_bw_sla 0.001 0.004 0.000 -0.003 0.000
AOU_WOA_clim 0.000 0.000 0.000 0.000 0.000
density_WOA_clim 0.000 0.000 0.026 -0.038 0.051
o2sat_WOA_clim 0.010 0.017 -0.004 -0.022 0.000
oxygen_WOA_clim 0.000 0.000 0.000 0.000 0.000
salinity_WOA_clim 0.000 0.000 0.000 0.000 0.000
conductivity_WOA_clim 0.000 0.000 0.000 0.000 0.000
nitrate_WOA_clim 0.000 0.000 0.000 0.000 0.000
phosphate_WOA_clim 0.000 0.000 0.052 -0.132 0.012
silicate_WOA_clim -0.016 0.000 -0.022 -0.001 -0.045
par 0.002 0.004 0.000 0.002 -0.002

Figure 11

Figure 11 can be produced using this code:

mats = list(alpha = cvres$bestres$alpha, beta = cvres$bestres$beta)

## Plot covariates
allnames = unique(c(colnames(mats$alpha), rownames(mats$beta[[1]])))
allnames = allnames[-which(allnames == "intp")]
Xsmall = cvres$bestres$X %>% as_tibble %>% select(one_of(allnames)) %>% select(sort(names(.)))
inds = round(seq(from = 1, to = ncol(Xsmall), length=5))

##' Make ticks from rownames of res$X (only for 1d data).
##' @param res Object of class |flowmix|.
add_date_ticks <- function(res){
  dates = sapply(as.Date(rownames(res$X)), 
                 toString)

  nums = as.numeric(as.factor(dates))
  left_ticks = sapply(sort(unique(nums)),function(ii){min(which(nums==ii))})
  left_ticks = c(left_ticks, res$TT)
  mid_ticks = sapply(sort(unique(nums)),function(ii){mean(which(nums==ii))})
  dates_mid_ticks = dates[round(mid_ticks)]
  axis(1, at=left_ticks, labels=FALSE)
  axis(1, at=mid_ticks, labels = dates_mid_ticks, tick=FALSE, las=2)
  axis(2)
}


## pdf(file = file.path(figdir, "1d-covariates.pdf"), width=20, height=15)
par(mfrow=c(2,2))
for(ii in 1:4){
  ind = (inds[ii]+1):inds[ii+1]
  cols = RColorBrewer::brewer.pal(length(ind), "Set2")
  matplot(Xsmall, axes=FALSE, col='grey', type='l', lwd=.5, lty=1, ylab="",xlab="")
  add_date_ticks(cvres$bestres)
  matlines(Xsmall[,ind], col=cols, lwd=3, lty=1)
  legend("topright", col=cols, lwd=3, lty=1, legend = colnames(Xsmall)[ind])
}

## graphics.off()

Figure 13

Figure 13 is produced using the second block (starting with # A 5 x 5 version of this):

load(file.path(cvres$destin, "meta.Rdata")) ## load ylist
## png(file=file.path(figdir, "1d-all-models.png"), width=5000, height=4000)
cv_gridsize = 10
par(mfrow = c(cv_gridsize, cv_gridsize))
for(ialpha in 1:cv_gridsize){
  for(ibeta in 1:cv_gridsize){
    res = cvres$bestreslist[[paste0(ialpha, "-", ibeta)]]
    res$prob = res$pie ## For back-compatibility; prob was called pie
    plot_1d(ylist = ylist, countslist = countslist, res = res,
            cex_data = 1)
    title(main = c(paste0("in-sample: ", round(min(res$objective),3)), "   ",
                   paste0("CV score : ", round(cvres$cvscore.mat[ialpha,ibeta],3))),
          cex.main=3)
    if(all(c(ialpha, ibeta) == cvres$min.ind)) box(lwd=10,col='blue')
  }
}
## graphics.off()

## A 5 x 5 version of this.
load(file.path(cvres$destin, "meta.Rdata")) ## load ylist
## png(file=file.path(figdir, "1d-all-models-reduced.png"), width=5000/2, height=3000/2)
cv_gridsize = 10
par(mfrow = c(cv_gridsize/2, cv_gridsize/2 ))
par(mar=c(0,0,5,0))
for(ialpha in c(1,3,5,7,10)){
  for(ibeta in c(1,4,6,8,10)){
    print(c(ialpha, ibeta))
    res = cvres$bestreslist[[paste0(ialpha, "-", ibeta)]]
    res$prob = res$pie ## For back-compatibility; prob was called pie
    plot_1d(ylist = ylist, countslist = countslist, res = res,
            no_axis=TRUE, cex_clust_label=NULL, omit_label=TRUE,
            omit_band=TRUE,
            cex_data = 1.7)
    title(main = c(paste0("In-sample: ", round(min(res$objective),3)), "   ",
                   paste0("CV score : ", round(cvres$cvscore.mat[ialpha,ibeta],3))),
          cex.main=4)##3.5)
    if(all(c(ialpha, ibeta) == cvres$min.ind)) box(lwd=10,col='blue')
  }
}
## graphics.off()

3d 10-cluster model results

### Load the model 
cvres_filename = paste0("cvres-", paste0(2,"-", 64,"-", 10), ".Rdata")

### cvres = blockcv_summary(2, 64, 10, 7, nrep = 5) 

### save(cvres, file = file.path(datadir, cvres_filename)) 

load(file = file.path(datadir, cvres_filename), verbose = FALSE)
res = cvres$bestres
destin = cvres$destin

## Quick adjustments, for back-compatibility
class(res) = "flowmix" ## for back-compatibility
res$prob = res$pie ## for back-compatibility

Table 3,4,5,6

Table 3,4,5,6 are produced using this code:

for(iclust in 1:10){
  onebeta = res$beta[[iclust]]
  colnames(onebeta) = c("diam", "red", "orange")
  knitr::kable(onebeta, digits=10, caption = paste0("Beta coefficientns for Cluster ", iclust)) %>% print()
}
Beta coefficientns for Cluster 1
diam red orange
intp 5.968095531 4.259373072 7.312468609
b1 0.000000000 0.000000000 0.000000000
b2 0.000000000 0.000000000 0.000000000
p1 0.000000000 0.088480468 0.000000000
p2 0.003237397 0.000000000 0.000000000
p3 0.002195150 0.004950264 -0.005398967
p4 0.000000000 -0.012164948 0.023084113
sss 0.000000000 0.000000000 0.000000000
sst 0.000000000 0.000000000 0.000000000
Fe 0.000000000 0.000000000 0.000000000
PP 0.000000000 0.000000000 0.000000000
Si 0.000000000 0.000000000 0.000000000
NO3 0.000000000 0.000000000 0.000000000
CHL 0.000000000 0.000000000 0.000000000
PHYC 0.000000000 0.000000000 0.000000000
PO4 0.000000000 0.000000000 0.000000000
O2 0.000000000 0.000000000 0.000000000
vgosa 0.000000000 0.000000000 0.000000000
vgos 0.000000000 0.000000000 0.000000000
sla 0.000000000 -0.001524180 0.025666353
ugosa 0.000000000 -0.044216683 0.000000000
ugos 0.000000000 0.000000000 0.000000000
wind_stress 0.000000000 0.000000000 0.000000000
eastward_wind 0.000000000 0.000000000 0.000000000
surface_downward_eastward_stress 0.000000000 0.000000000 0.000000000
wind_speed 0.000000000 0.000000000 0.000000000
surface_downward_northward_stress 0.000000000 0.000000000 0.000000000
northward_wind 0.000000000 0.000000000 0.000000000
ftle_bw_sla 0.000000000 0.000000000 0.000000000
disp_bw_sla 0.000000000 0.000000000 0.000000000
AOU_WOA_clim 0.000000000 0.000000000 0.000000000
density_WOA_clim 0.000000000 0.000000000 0.000000000
o2sat_WOA_clim 0.000000000 0.000000000 0.000000000
oxygen_WOA_clim 0.000000000 0.000000000 0.000000000
salinity_WOA_clim 0.000000000 0.000000000 0.000000000
conductivity_WOA_clim 0.000000000 0.000000000 0.000000000
nitrate_WOA_clim 0.000000000 0.000000000 0.000000000
phosphate_WOA_clim 0.000000000 0.000000000 0.000000000
silicate_WOA_clim 0.000000000 0.000000000 0.000000000
par 0.000000000 0.000000000 0.000000000
Beta coefficientns for Cluster 2
diam red orange
intp 5.9284345966 4.85418418 2.81943125
b1 0.0000000000 0.00000000 0.00000000
b2 0.0000000000 0.00000000 0.00000000
p1 -0.0001051634 0.00000000 0.00000000
p2 0.0000000000 0.00000000 0.00000000
p3 0.0000000000 0.00000000 0.00000000
p4 0.0000000000 0.00000000 0.00000000
sss 0.0000000000 0.00000000 0.00000000
sst 0.0000000000 0.00000000 0.00000000
Fe 0.0000000000 0.00000000 0.07211896
PP 0.0000000000 0.00000000 0.00000000
Si 0.0000000000 0.00000000 0.00000000
NO3 0.0000000000 0.00000000 0.00000000
CHL 0.0000000000 0.00000000 0.00000000
PHYC 0.0000000000 0.01947219 0.00000000
PO4 0.0000000000 0.03782072 0.00000000
O2 0.0000000000 0.00000000 0.00000000
vgosa 0.0000000000 0.00000000 0.00000000
vgos 0.0000000000 0.00000000 0.00000000
sla 0.0000000000 0.00000000 0.04791567
ugosa 0.0000000000 0.00000000 0.00000000
ugos 0.0000000000 0.00000000 0.00000000
wind_stress 0.0000000000 0.00000000 0.00000000
eastward_wind 0.0000000000 0.00000000 0.00000000
surface_downward_eastward_stress 0.0000000000 0.00000000 0.00000000
wind_speed 0.0000000000 0.00000000 0.00000000
surface_downward_northward_stress 0.0000000000 0.00000000 0.00000000
northward_wind 0.0000000000 0.00000000 0.00000000
ftle_bw_sla 0.0000000000 0.00000000 0.00000000
disp_bw_sla 0.0000000000 0.00000000 0.00000000
AOU_WOA_clim 0.0000000000 0.00000000 0.00000000
density_WOA_clim 0.0000000000 0.00000000 0.00000000
o2sat_WOA_clim -0.0051724217 0.00000000 0.00000000
oxygen_WOA_clim 0.0000000000 0.00000000 0.00000000
salinity_WOA_clim 0.0000000000 0.00000000 0.00000000
conductivity_WOA_clim 0.0000000000 0.00000000 0.00000000
nitrate_WOA_clim 0.0000000000 0.00000000 0.00000000
phosphate_WOA_clim 0.0000000000 0.00000000 0.00000000
silicate_WOA_clim 0.0000000000 0.00000000 -0.01639025
par 0.0000000000 0.00000000 0.00000000
Beta coefficientns for Cluster 3
diam red orange
intp 3.053828688 2.734432867 3.066667492
b1 0.000000000 0.000000000 0.188918929
b2 -0.030552519 -0.214914300 -0.046613588
p1 0.028770267 -0.016449035 0.000000000
p2 0.033964054 -0.005357155 0.017515021
p3 0.039639725 0.022252885 0.000000000
p4 0.000000000 0.018108747 0.024554493
sss 0.000000000 0.026642658 0.000000000
sst 0.000000000 0.000000000 -0.001066180
Fe -0.028930710 0.000000000 0.000000000
PP 0.000000000 -0.113621472 -0.099036941
Si 0.266369620 0.000000000 0.000000000
NO3 0.000000000 0.000000000 0.000000000
CHL -0.137182964 0.184773277 0.081165090
PHYC 0.000000000 0.000000000 0.000000000
PO4 0.000000000 0.000000000 0.000000000
O2 0.000000000 0.000000000 0.034039716
vgosa 0.000000000 -0.019974801 0.000000000
vgos 0.043537783 0.000000000 0.000000000
sla 0.049054132 -0.011223886 0.014263833
ugosa 0.000000000 0.000000000 0.000000000
ugos -0.073874323 0.000000000 0.062561257
wind_stress 0.000000000 -0.096102928 -0.056282580
eastward_wind -0.016324939 0.000000000 0.000000000
surface_downward_eastward_stress 0.000000000 -0.050850935 0.000000000
wind_speed -0.111971670 0.077574250 0.097343356
surface_downward_northward_stress 0.021993977 0.000000000 0.000000000
northward_wind 0.000000000 0.014900891 -0.001973156
ftle_bw_sla 0.098820130 0.075252460 0.006401732
disp_bw_sla 0.066465241 0.000000000 -0.015352797
AOU_WOA_clim 0.000000000 0.000000000 -0.168113341
density_WOA_clim 0.000000000 0.000000000 0.019794785
o2sat_WOA_clim 0.001731463 0.135285232 0.000000000
oxygen_WOA_clim 0.000000000 0.000000000 0.000000000
salinity_WOA_clim 0.000000000 0.000000000 0.000000000
conductivity_WOA_clim 0.000000000 0.000000000 0.000000000
nitrate_WOA_clim 0.000000000 -0.017238285 0.000000000
phosphate_WOA_clim 0.000000000 0.000000000 0.042338449
silicate_WOA_clim 0.000000000 0.000000000 -0.100771898
par 0.019632919 0.000000000 0.039190904
Beta coefficientns for Cluster 4
diam red orange
intp 3.94478275 0.800939334 0.5510790298
b1 0.00000000 0.000000000 0.0000000000
b2 0.00000000 -0.011203322 0.0470913431
p1 0.00000000 0.000000000 -0.0051267445
p2 0.00000000 0.000000000 -0.0003768198
p3 0.00000000 0.000000000 -0.0046129231
p4 0.00000000 0.000000000 -0.0049090076
sss 0.00000000 0.000000000 -0.0199267703
sst 0.00000000 0.000000000 0.0000000000
Fe 0.00000000 0.000000000 0.0000000000
PP 0.00000000 0.000000000 0.0072871711
Si 0.00000000 0.000000000 0.0000000000
NO3 0.00000000 0.000000000 0.0000000000
CHL 0.00000000 0.000000000 0.0000000000
PHYC 0.00000000 0.000000000 0.0000000000
PO4 0.00000000 0.000000000 0.0015115451
O2 0.00000000 0.000000000 0.0000000000
vgosa 0.00000000 0.000000000 0.0000000000
vgos 0.00000000 0.000000000 0.0000000000
sla 0.00000000 0.000000000 0.0000000000
ugosa 0.00000000 -0.005928616 0.0000000000
ugos 0.00000000 0.000000000 0.0000000000
wind_stress 0.00000000 0.000000000 0.0000000000
eastward_wind 0.00000000 0.044257257 0.0000000000
surface_downward_eastward_stress 0.00000000 0.000000000 0.0156934800
wind_speed 0.08987578 0.000000000 0.0000000000
surface_downward_northward_stress 0.00000000 0.000000000 0.0000000000
northward_wind 0.00000000 -0.001137329 0.0000000000
ftle_bw_sla 0.00000000 0.000000000 0.0000000000
disp_bw_sla 0.00000000 0.000000000 0.0000000000
AOU_WOA_clim 0.00000000 0.000000000 0.0000000000
density_WOA_clim 0.00000000 0.000000000 0.0000000000
o2sat_WOA_clim 0.00000000 0.000000000 0.0000000000
oxygen_WOA_clim 0.00000000 0.000000000 0.0000000000
salinity_WOA_clim 0.00000000 0.000000000 0.0000000000
conductivity_WOA_clim 0.00000000 0.000000000 0.0000000000
nitrate_WOA_clim 0.00000000 0.000000000 0.0445281005
phosphate_WOA_clim 0.01330943 0.000000000 0.0000000000
silicate_WOA_clim 0.00000000 0.000000000 0.0000000000
par 0.00000000 0.000000000 -0.0022121860
Beta coefficientns for Cluster 5
diam red orange
intp 4.986441529 5.324611782 0.4936000916
b1 0.000000000 0.000000000 0.0056939175
b2 0.062266024 -0.156980029 0.0000000000
p1 0.000000000 -0.003904175 0.0000000000
p2 0.041519223 -0.001549597 -0.0020668092
p3 0.035577815 0.034980200 -0.0062915949
p4 0.000000000 0.065074590 -0.0027824560
sss 0.000000000 0.000000000 0.0000000000
sst 0.000000000 0.000000000 0.0000000000
Fe 0.000000000 0.000000000 0.0000000000
PP 0.000000000 -0.395517363 0.0000000000
Si 0.000000000 0.000000000 0.0000000000
NO3 0.000000000 0.000000000 0.0000000000
CHL 0.000000000 0.214728693 0.0389324522
PHYC 0.000000000 0.000000000 0.0000000000
PO4 0.000000000 0.000000000 0.0000000000
O2 0.000000000 0.000000000 0.0000000000
vgosa 0.093737891 0.000000000 0.0108960904
vgos 0.000000000 0.000000000 0.0000000000
sla 0.000000000 0.000000000 0.0000000000
ugosa 0.039513833 0.000000000 0.0048786632
ugos 0.000000000 0.013910597 0.0000000000
wind_stress 0.000000000 -0.062007723 0.0000000000
eastward_wind 0.000000000 0.000000000 0.0184339923
surface_downward_eastward_stress 0.000000000 0.000000000 0.0000000000
wind_speed -0.092006661 0.000000000 0.0099805389
surface_downward_northward_stress 0.002308997 -0.072639723 0.0000000000
northward_wind 0.000000000 0.000000000 0.0000000000
ftle_bw_sla 0.000000000 0.056271420 0.0000000000
disp_bw_sla 0.000000000 0.000000000 0.0000000000
AOU_WOA_clim 0.000000000 0.000000000 0.0000000000
density_WOA_clim 0.046191169 0.000000000 0.0000000000
o2sat_WOA_clim -0.175196083 0.000000000 -0.0280276234
oxygen_WOA_clim 0.000000000 0.000000000 0.0000000000
salinity_WOA_clim 0.000000000 0.000000000 0.0000000000
conductivity_WOA_clim 0.000000000 0.000000000 0.0000000000
nitrate_WOA_clim 0.000000000 0.048632912 0.0000000000
phosphate_WOA_clim 0.128105135 0.096144750 0.0275039300
silicate_WOA_clim 0.000000000 0.095391721 -0.0019678715
par 0.034601411 -0.016013761 0.0007855167
Beta coefficientns for Cluster 6
diam red orange
intp 3.482483385 1.5918171 1.532339
b1 0.000000000 0.0000000 0.000000
b2 0.000000000 0.0000000 0.000000
p1 0.000000000 0.0000000 0.000000
p2 -0.028452219 0.0000000 0.000000
p3 0.000000000 0.0000000 0.000000
p4 0.000000000 0.0000000 0.000000
sss 0.000000000 0.0000000 0.000000
sst 0.000000000 0.0000000 0.000000
Fe 0.000000000 0.0000000 0.000000
PP 0.000000000 0.0000000 0.000000
Si 0.000000000 0.0000000 0.000000
NO3 0.000000000 0.0000000 0.000000
CHL 0.000000000 0.0000000 0.000000
PHYC 0.000000000 0.0000000 0.000000
PO4 0.000000000 0.0000000 0.000000
O2 0.000000000 0.0000000 0.000000
vgosa 0.071386211 0.0000000 0.000000
vgos 0.000000000 0.0000000 0.000000
sla 0.002178113 0.0000000 0.000000
ugosa 0.000000000 0.0000000 0.000000
ugos 0.000000000 0.0000000 0.000000
wind_stress 0.000000000 0.0000000 0.000000
eastward_wind 0.000000000 0.0443385 0.000000
surface_downward_eastward_stress 0.000000000 0.0000000 0.000000
wind_speed 0.000000000 0.0000000 0.000000
surface_downward_northward_stress 0.000000000 0.0000000 0.000000
northward_wind 0.025333431 0.0000000 0.000000
ftle_bw_sla 0.000000000 0.0000000 0.000000
disp_bw_sla 0.000000000 0.0000000 0.000000
AOU_WOA_clim 0.000000000 0.0000000 0.000000
density_WOA_clim 0.000000000 0.0000000 0.000000
o2sat_WOA_clim 0.000000000 0.0000000 0.000000
oxygen_WOA_clim 0.000000000 0.0000000 0.000000
salinity_WOA_clim 0.000000000 0.0000000 0.000000
conductivity_WOA_clim 0.000000000 0.0000000 0.000000
nitrate_WOA_clim 0.004990995 0.0000000 0.000000
phosphate_WOA_clim 0.197216611 0.0000000 0.000000
silicate_WOA_clim 0.000000000 0.0000000 0.000000
par 0.000000000 0.0000000 0.000000
Beta coefficientns for Cluster 7
diam red orange
intp 6.468244652 6.690494715 1.457910313
b1 0.000000000 -0.010066972 0.000000000
b2 0.164754781 0.000000000 0.039131959
p1 0.031037845 0.000000000 0.000000000
p2 0.043703004 -0.003875314 -0.005956309
p3 0.042477276 0.000000000 -0.000941041
p4 0.009046667 0.000000000 -0.011118713
sss -0.049222368 0.000000000 0.000000000
sst 0.000000000 0.000000000 0.000000000
Fe 0.000000000 0.000000000 0.000000000
PP 0.269983557 -0.052958459 0.022525970
Si 0.000000000 0.030137957 0.000000000
NO3 0.000000000 0.000000000 0.000000000
CHL -0.392623178 0.000000000 0.000000000
PHYC 0.000000000 0.000000000 0.000000000
PO4 0.000000000 0.000000000 0.000000000
O2 0.000000000 0.000000000 0.000000000
vgosa 0.003715129 0.000000000 0.000000000
vgos 0.059892049 0.000000000 0.000000000
sla 0.000000000 -0.039857430 0.000000000
ugosa 0.000000000 0.000000000 0.000000000
ugos 0.014007039 0.000000000 0.000000000
wind_stress 0.000000000 0.002194998 0.004492496
eastward_wind 0.000000000 0.000000000 0.067508265
surface_downward_eastward_stress 0.000000000 0.014085674 0.000000000
wind_speed -0.021648796 0.024419466 0.000000000
surface_downward_northward_stress 0.000000000 0.000000000 0.000000000
northward_wind 0.030869437 -0.018596375 0.000000000
ftle_bw_sla 0.004383963 0.000000000 0.000000000
disp_bw_sla 0.006809223 0.000000000 0.000000000
AOU_WOA_clim 0.000000000 0.000000000 0.000000000
density_WOA_clim 0.000000000 0.000000000 0.000000000
o2sat_WOA_clim 0.118108706 0.136655866 0.000000000
oxygen_WOA_clim 0.000000000 0.000000000 0.000000000
salinity_WOA_clim 0.000000000 0.000000000 0.000000000
conductivity_WOA_clim 0.000000000 0.000000000 0.000000000
nitrate_WOA_clim -0.008558563 0.000000000 -0.051040858
phosphate_WOA_clim 0.000000000 0.000000000 0.000000000
silicate_WOA_clim -0.032048750 0.008768610 -0.017496671
par 0.040306274 -0.015184376 0.056219340
Beta coefficientns for Cluster 8
diam red orange
intp 2.067859574 0.6553069324 0.3586214454
b1 0.000000000 -0.0303536773 0.0000000000
b2 0.030442408 -0.0485203544 0.0000000000
p1 -0.016742573 0.0121478092 -0.0015608202
p2 -0.007582239 0.0001740984 -0.0012395712
p3 -0.026307710 0.0000000000 -0.0030449855
p4 -0.006150394 0.0000000000 -0.0004985803
sss 0.000000000 0.0000000000 0.0000000000
sst 0.000000000 -0.0359499740 0.0000000000
Fe 0.000000000 0.0000000000 0.0000000000
PP 0.000000000 -0.0156779941 0.0000000000
Si -0.009905474 0.0000000000 -0.0015481463
NO3 0.000000000 0.0000000000 0.0000000000
CHL 0.000000000 0.0000000000 0.0005304632
PHYC 0.033957894 0.0000000000 0.0000000000
PO4 0.000000000 0.0000000000 0.0000000000
O2 0.000000000 0.0000000000 0.0187033314
vgosa 0.060682121 0.0000000000 0.0000000000
vgos 0.000000000 0.0019360080 0.0020665502
sla 0.073511978 -0.0022250961 0.0000000000
ugosa 0.000000000 0.0000000000 0.0030580396
ugos 0.034208312 0.0000000000 0.0000000000
wind_stress 0.000000000 -0.0319185848 0.0000000000
eastward_wind 0.113847926 0.0098551581 0.0113342962
surface_downward_eastward_stress 0.138645188 0.0000000000 0.0000000000
wind_speed 0.135558314 0.0432119919 0.0052366460
surface_downward_northward_stress 0.005789204 0.0000000000 0.0039504788
northward_wind 0.000000000 -0.0120681600 0.0000000000
ftle_bw_sla 0.000000000 0.0000000000 0.0015998364
disp_bw_sla -0.056236936 0.0000000000 0.0000000000
AOU_WOA_clim 0.000000000 0.0000000000 0.0000000000
density_WOA_clim 0.030581519 0.0000000000 0.0000000000
o2sat_WOA_clim 0.038474202 0.0250866785 -0.0004554307
oxygen_WOA_clim 0.000000000 0.0000000000 0.0000000000
salinity_WOA_clim 0.000000000 0.0000000000 0.0000000000
conductivity_WOA_clim 0.000000000 0.0000000000 0.0000000000
nitrate_WOA_clim 0.000000000 0.0000000000 0.0000000000
phosphate_WOA_clim 0.000000000 0.0000000000 0.0194514971
silicate_WOA_clim -0.121452678 -0.0026934780 -0.0048742566
par 0.000000000 0.0009309060 -0.0003262723
Beta coefficientns for Cluster 9
diam red orange
intp 6.890361668 1.52187313 1.743250242
b1 0.000000000 0.00000000 0.000000000
b2 0.000000000 0.00000000 0.000000000
p1 0.000000000 0.00000000 0.000000000
p2 0.000000000 0.00000000 0.000000000
p3 0.000000000 0.00000000 0.000000000
p4 0.000000000 0.00000000 0.000000000
sss 0.000000000 0.00000000 -0.006719930
sst 0.011253357 0.00000000 0.000000000
Fe 0.000000000 -0.01881572 0.000000000
PP 0.000000000 0.00000000 0.077935649
Si 0.000000000 0.00000000 0.000000000
NO3 0.000000000 0.00000000 0.000000000
CHL 0.000000000 0.00000000 0.000000000
PHYC 0.000000000 0.00000000 0.000000000
PO4 0.000000000 0.00000000 0.000000000
O2 0.000000000 0.00000000 0.000000000
vgosa 0.000000000 0.00000000 0.000000000
vgos 0.000000000 0.00000000 0.000000000
sla 0.000000000 0.00000000 0.000000000
ugosa 0.000000000 0.00000000 0.000000000
ugos 0.000000000 0.00000000 0.000000000
wind_stress 0.000000000 0.00000000 0.000000000
eastward_wind -0.009923210 0.00000000 0.000000000
surface_downward_eastward_stress 0.000000000 0.00000000 0.000000000
wind_speed 0.000000000 0.00000000 0.000000000
surface_downward_northward_stress 0.000000000 0.00000000 0.000000000
northward_wind 0.000000000 0.00000000 0.000000000
ftle_bw_sla 0.000000000 0.00000000 0.000000000
disp_bw_sla 0.003051813 -0.01542399 0.000000000
AOU_WOA_clim 0.000000000 0.00000000 0.000000000
density_WOA_clim 0.000000000 0.00000000 0.000000000
o2sat_WOA_clim -0.046664954 0.01486046 0.000000000
oxygen_WOA_clim 0.000000000 0.00000000 0.000000000
salinity_WOA_clim 0.000000000 0.00000000 0.000000000
conductivity_WOA_clim 0.000000000 0.00000000 0.000000000
nitrate_WOA_clim 0.000000000 0.00000000 0.000000000
phosphate_WOA_clim 0.000000000 0.00000000 0.000000000
silicate_WOA_clim 0.000000000 0.02343054 0.000000000
par 0.000000000 -0.01471873 -0.006898926
Beta coefficientns for Cluster 10
diam red orange
intp 1.415990835 1.362324973 0.3568267870
b1 0.171902625 0.000000000 0.0000000000
b2 0.049009117 -0.129168463 -0.0055911873
p1 0.040388261 -0.019673175 0.0008794797
p2 0.072969681 0.012931907 0.0000000000
p3 0.022498169 -0.027850210 -0.0000365907
p4 -0.002910622 0.000000000 0.0000000000
sss -0.054764821 0.000000000 0.0000000000
sst -0.314879406 0.000000000 0.0000000000
Fe 0.000000000 -0.006268936 0.0000000000
PP 0.076559334 0.000000000 0.0000000000
Si 0.000000000 0.014480287 0.0000000000
NO3 0.000000000 0.000000000 0.0000000000
CHL 0.000000000 0.098410546 0.0000000000
PHYC 0.006361228 0.000000000 0.0000000000
PO4 0.000000000 0.000000000 0.0000000000
O2 0.000000000 0.000000000 0.0000000000
vgosa 0.000000000 -0.024449143 0.0000000000
vgos -0.020260243 0.000000000 -0.0003841106
sla -0.027651693 -0.013693736 0.0000000000
ugosa -0.046100353 -0.015591432 0.0003389835
ugos 0.056819468 0.000000000 0.0000000000
wind_stress -0.018279871 0.000000000 0.0000000000
eastward_wind 0.000000000 0.047565529 0.0122879348
surface_downward_eastward_stress -0.028759822 -0.083413614 0.0000000000
wind_speed 0.000000000 0.003169275 0.0072140418
surface_downward_northward_stress 0.000000000 0.000000000 0.0000000000
northward_wind 0.001262733 0.000000000 0.0000000000
ftle_bw_sla 0.002712154 -0.020627653 0.0000000000
disp_bw_sla -0.017695401 -0.007581416 -0.0006605134
AOU_WOA_clim 0.000000000 -0.029385509 0.0000000000
density_WOA_clim -0.060388349 0.019958541 0.0000000000
o2sat_WOA_clim -0.040655042 0.000000000 0.0000000000
oxygen_WOA_clim 0.000000000 0.000000000 0.0000000000
salinity_WOA_clim 0.000000000 0.000000000 -0.0087392612
conductivity_WOA_clim 0.000000000 0.000000000 0.0000000000
nitrate_WOA_clim -0.181959755 -0.223673418 0.0000000000
phosphate_WOA_clim -0.086413882 0.000000000 0.0000000000
silicate_WOA_clim -0.040446198 -0.003833241 0.0000000000
par 0.024380325 -0.059013480 0.0000000000 
t(res$alpha) %>%  knitr::kable(digits=10, caption = paste0("Alpha coefficients")) %>% print()
Alpha coefficients
clust-1 clust-2 clust-3 clust-4 clust-5 clust-6 clust-7 clust-8 clust-9 clust-10
intp -5.354472250 -3.13822026 -2.48933025 -2.79681799 -2.37744920 -3.312255 -1.85996671 -1.913991746 -3.07497797 -2.30878410
b1 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
b2 0.000000000 0.00000000 0.00000000 0.00000000 -0.01726096 0.000000 0.00000000 0.716205823 0.00000000 0.00000000
p1 0.085836335 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.011783557 0.00000000 0.00000000
p2 0.000000000 0.00000000 0.00000000 0.00000000 -0.12997729 0.000000 0.00000000 0.000000000 0.00000000 0.05748439
p3 0.000000000 0.00000000 0.04719697 0.00000000 0.00000000 0.000000 0.00000000 -0.038274581 0.00000000 0.01806784
p4 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 -0.023251193 0.00000000 0.00000000
sss 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.02290001
sst 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.86791695
Fe 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
PP 0.000000000 0.00000000 0.19469763 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
Si 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
NO3 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
CHL 0.000000000 0.00000000 0.00000000 0.00000000 0.21432808 0.000000 0.00000000 0.000000000 0.00000000 -0.06066823
PHYC 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
PO4 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
O2 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
vgosa 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 -0.060801105 0.00000000 0.00000000
vgos -0.002114332 0.00000000 -0.08849818 0.00000000 0.00000000 0.000000 0.07933362 0.000000000 0.00000000 0.00000000
sla 0.000000000 0.02483285 -0.05579800 0.00000000 0.20098887 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
ugosa 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.02540237 0.000000000 0.00000000 -0.03350112
ugos 0.000000000 0.00000000 0.17528502 0.00000000 -0.18775088 0.000000 0.00000000 0.000000000 0.00000000 -0.06048606
wind_stress -0.436012754 0.03584415 -0.37055443 0.00000000 -0.08539663 0.000000 0.00000000 0.107894170 0.16017099 0.00000000
eastward_wind 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
surface_downward_eastward_stress 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.140214981 0.00000000 0.00000000
wind_speed 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 -0.25457259
surface_downward_northward_stress 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
northward_wind 0.000000000 0.00000000 -0.18059294 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.01933629 -0.10500179
ftle_bw_sla 0.000000000 0.00000000 -0.18030660 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
disp_bw_sla 0.000000000 0.00000000 0.00000000 0.00000000 -0.16096143 0.000000 0.00000000 0.000000000 0.00000000 0.01098454
AOU_WOA_clim 0.000000000 0.00000000 0.00000000 0.03189847 -0.01966688 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
density_WOA_clim 0.000000000 0.00000000 0.00000000 0.00000000 0.42886047 0.000000 0.11105708 0.000000000 0.00000000 0.00000000
o2sat_WOA_clim 0.000000000 0.00000000 0.45542258 -0.04876088 0.00000000 0.000000 0.00000000 0.002985964 -0.14259323 0.00000000
oxygen_WOA_clim 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
salinity_WOA_clim 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
conductivity_WOA_clim 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
nitrate_WOA_clim 0.000000000 0.29802627 -0.34708008 0.00000000 0.00000000 0.000000 0.00000000 -0.005186623 0.01487063 0.00000000
phosphate_WOA_clim 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.34796662 0.000000000 0.00000000 -0.93665103
silicate_WOA_clim 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.000000000 0.00000000 0.00000000
par 0.431517342 0.00000000 0.00000000 0.00000000 0.00000000 0.000000 0.00000000 0.012834378 0.00000000 -0.07906631

Figure 8

Figure 8 is produced using this code:

## Rescale diameters (not used anymore since diam in the data is in original
## range, and is shifted/scaled to be in the same scale as CHL and PE only at 
## simulation time.)
source("gradients2-helpers.R")

## Temporary *undoing* of the scaling in the estimated model.
range.diam = c(-0.6706111,  1.3070264)
range.chl =  c(0.101313, 7.955767)
range.pe = c(0.101313, 7.955767)
  
## Modify (rescale/shift) the res$mn
width.diam = range.diam[2] - range.diam[1]
idim = 1
res$mn = res$mn[,idim,, drop=FALSE] / max(range.chl) * width.diam + min(range.diam)

## Undoing the scaling in res$sigma
res$sigma = cvres$bestres$sigma[,idim, idim, drop=FALSE] 
myscale = max(range.chl) * width.diam
res$sigma = res$sigma / myscale

## Load the data
datobj = readRDS(file.path(datadir, "MGL1704-hourly-paper-1d-diam.RDS"))
ylist = datobj$ylist
countslist = datobj$countslist

## Make the 1d data+model plot
res$dimdat = 1
cols = plot_1d(ylist = ylist, countslist = countslist, res, scale=TRUE, omit_band=TRUE, omit_label=TRUE)

## Additionally, add y axis label
mtext("Diameter", side=2, line=2.5)

## Additionally, write the cluster numbers with the same colors
res <- reorder_clust(res)
TT = length(ylist)
for(iclust in 1:res$numclust){
  y = res$mn[1, 1, iclust]
  if(iclust==3) y = y + 0.05
  if(iclust==4) y = y - 0.05
  text(x = res$TT/50 - 5, y = y,
       label = paste0("Cluster ", iclust), cex = 1, col="black")##col=cols[iclust])
}

It’s worth explaining a bit about data rescaling for Gradients 2. The three dimensions of the original cytogram data (e.g. MGL1704-hourly-paper.Rdata or in MGL1704-hourly.Rdata) all have a range roughly between 0 and 8, as you can see here:

datobj = readRDS(file = file.path(datadir, "MGL1704-hourly-paper.RDS"))
datobj$ylist[[2]] %>% apply(2, range) %>% knitr::kable() %>% print()
diam_mid chl_small pe
0.000000 0.101313 0.101313
7.955767 7.955767 7.955767

This is the result of shifting and stretching the diameter dimension (the first of the three dimensions), which originally had a range of:

datobj = readRDS(file.path(datadir, "MGL1704-hourly-paper-1d-diam.RDS"));
datobj$ylist[[1]] %>% range() %>% cbind() %>% `colnames<-`(c("diam_mid")) %>% knitr::kable() %>% print()
diam_mid
-0.6706111
1.3070264

to be the same range as the other two dimensions (Chl and PE), between \(0\) and \(7.955767\). So the code above is just reversing this shifting and scaling.

Figure 9

Figure 9 is produced using this code:

plot_3d_fourpanels <- function(obj, ## Understandably, data (ylist) might not be in the object.
                               ylist, countslist = NULL, ## The time point of interest, out of 1:TT
                               tt, ## Other options.
                               ## 2d scatterplot options
                               show.xb.constraint = FALSE, cex.fac.2d = 1, par_cex_2d = 1,
                               pt_col = rgb(0 ,0, 1, 0.1), ## 3d scatterplot options
                               cex.fac.3d = 1, ## 3d scatterplot options
                               destin = NULL){

  ## Define layout
  par(mfrow = c(1, 4))
  ## par(oma = c(3, 1, 2, 1)) ## setting outer margin

  ## Setup
  TT = length(ylist)
  assertthat::assert_that(tt %in% 1:TT)
  all.y = do.call(rbind, ylist)
  only_plot_cytograms = (is.null(obj))
  if(!only_plot_cytograms){
    obj = reorder_clust(obj)
    mns = obj$mn
    numclust = obj$numclust
    p = ncol(obj$X)
  }

  ## Scale the biomass (|countslist|) by the total biomass in that cytogram.
  counts_sum = sapply(countslist, sum)
  fac = median(counts_sum)
  countslist = lapply(countslist, function(counts)counts/sum(counts) * fac)


  ###############################
  ## Make the three data plots ##
  ###############################
  ## par(mar=c(1,1,3,1))
  par(mar = c(5.1, 5.1, 4.1, 2.1))
  dimslist = list(1:2, 2:3, c(3,1))
  for(dims in dimslist){
    ## one_dim_scatterplot(ylist, obj, tt,
    scatterplot_2d(ylist, obj, tt,
                   countslist = countslist,
                   dims = dims,
                   cex_fac = cex.fac.2d,
                   pt_col = pt_col,
                   lwd = 2)
  }
  par(cex=0.8)

  par(mar=c(1,1,3,1))
  ## phis = c(10,50)
  phis = c(10)
  for(phi in phis){
    one_3d_plot(ylist, obj, tt, countslist = countslist, phi = phi,
                cex.fac = cex.fac.3d)
  }
}

## Load 3d model
cvres_filename = paste0("cvres-", paste0(2,"-", 64,"-", 10), ".Rdata")
load(file = file.path(datadir, cvres_filename), verbose = FALSE)
res = cvres$bestres
class(res) = "flowmix" ## for back-compatibility
res$prob = res$pie ## for back-compatibility

datobj = readRDS(file = file.path(datadir, "MGL1704-hourly-paper.RDS"))
ylist = datobj$ylist
countslist = datobj$countslist
tt = 10
plot_3d_fourpanels(res, ylist, countslist, tt = tt, destin = destin,
                   cex.fac = 10, pt_col = rgb(0,0,1,0.2), cex.fac.3d = 0.4)

Figure 10

Figure 10 is produced using this code:

(The file prochloro-gates.Rdata comes from paper-code/gating.Rmd).

## Load the results from traditional gating of Prochlorococcus
load(file = file.path(datadir, "prochloro-gates.Rdata"))

## 1. First, make 3-minute resolution prochlorococcus counts/biomass
prop_biomass = (prochloro_biomass_list %>% unlist) / (total_biomass_list %>% unlist)
names(prop_biomass) = names(prochloro_biomass_list)

prop_count = (prochloro_count_list %>% unlist) / (total_count_list %>% unlist)
names(prop_count) = names(prochloro_count_list)

## 2. Now, hourly aggregation of prochlorococcus counts/biomass.
times = names(prochloro_biomass_list)
obj = get_membership(times) ## borrowed from the "cruisedat" package
membership = obj$membership
time = obj$time

hourly_prochloro_biomass <- sapply(membership, function(memb){
  prochloro_biomass_list[memb] %>% unlist() %>% sum() })
hourly_total_biomass <- sapply(membership, function(memb){
  total_biomass_list[memb] %>% unlist() %>% sum()  })
names(hourly_prochloro_biomass) = names(hourly_total_biomass) = time
hourly_prop_biomass = hourly_prochloro_biomass / hourly_total_biomass
## Make plot 
## pdf(file = file.path(figdir, "3d-prochloro-prob.pdf"),
##     width = 8,height = 4)
par(mar = c(5.1, 4.1, 2.1, 2.1))
cols = RColorBrewer::brewer.pal(10, "Paired")
iclusts = c(10)
flowmix::plot_prob(res, iclusts, main = "", cols = cols)
lines(hourly_prop_biomass[-c(1:48)],
      lty = 1, type = 'l',lwd = 2,
      col='grey40')
legend("topleft",
       lwd = c(3, rep(3, length(iclusts))),
       col=c('grey40', cols[iclusts]),
       lty = 1,
       legend=c("Gated Prochlorococcus\nbiomass proportion", paste0("Cluster ", iclusts, " probability")),
       bg="white")

## graphics.off()

Figure 14

Figure 14 is produced using this code:

## ## Load 3d model
## cvres_filename = paste0("cvres-", paste0(2,"-", 64,"-", 10), ".Rdata")
## load(file = file.path(datadir, cvres_filename), verbose = FALSE)
## res = cvres$bestres
## class(res) = "flowmix" ## for back-compatibility
## res$prob = res$pie ## for back-compatibility

## ## Load 3d data
## datobj = readRDS(file = file.path(datadir, "MGL1704-hourly-paper.RDS"))
## ylist = datobj$ylist
## countslist = datobj$countslist

tt = 10 ## WHAT is this?
plot_3d(res, ylist, countslist, tt = tt, cex.fac = 10, pt_col = rgb(0,0,1,0.2),
        cex.fac.3d = 0.4,
        lat = datobj$lat,
        lon = datobj$lon)